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Abstract 

We present a generalization of the Li, Nunes and Vanderbilt density-matrix 



^ ■ method to the case of a non-orthogonal set of basis functions. A representation 

in , 

of the real-space density matrix is chosen in such a way that only the overlap 
f~^ ■ matrix, and not its inverse, appears in the energy functional. The generalized 

0^ . energy functional is shown to be variational with respect to the elements of 



the density matrix, which typically remains well localized. 
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Recently, the search for so-called "order- A^" methods (for which the computational effort 
scales only linearly with system size A^) has led to the development of a number of new real- 
space approaches to the solution of the electronic structure problenJiTQ. These are based 
either on the use of a localized, Wannier-like representation of the occupied subspaceH'S, or 
on the locality of the real-space density matrixp'a. In the latter case an energy functional is 
defined such that the variational degrees of freedom are the matrix elements of the density 
matrix in a real-space-localized set of orthonormal orbitals. However, in many situations 
it is more convenient to work with a non-orthogonal basis (e.g., LCAO calculations using 
Gaussian orbitals). For that reason, it becomes desirable to extend the density- matrix based 
approaches to those cases. 

In this paper we show how the approach proposed by Li, Nunes and Vanderbilt (LNV)a 
can be extended to a non-orthogonal basis. This is done by introducing a new quantity 
p = S^^pS^^ (where S is the overlap matrix Sij = {(pi \ (j)j), and p = {(pi \ p \ (pj) is the 
density matrix) as an alternative representation for the density operator, which is shown to 
have similar localization properties as p. Using p we write a generalized expression for the 
total energy in which the inverse overlap matrix S~^ does not appear explicitly; moreover 
the generalized-density-matrix (GDM) functional is shown to be variational with respect to 

P- 

First, we briefly review the LNV approach as applied to an orthonormal basisa. For 
simplicity we consider a tight-binding description of a system formed by replicating a unit 
(super)cell containing A^ atoms with M basis orbitals per site. For the moment we assume 
that the basis orbitals {(pi} are orthonormal, i.e., {(pi \ (pj) = 6ij. For the eigenstates of the 
Hamiltonian, H \ ipn) = Cn | i'n), we write 



V'n) =J2'^rii\ (pi) ■ (1) 



The density matrix is defined as 



Pij ~ Z^ '^ni '^nj ) \^) 



where i and j run over all tight-binding basis orbitals in the system and n runs over the 
occupied eigenstates of H. Recall that p is a projection onto the subspace of occupied states, 
and therefore it obeys the idempotency requirement p^ = p. 

As discussed in Ref.lil, both the standard k-space diagonalization of H and the minimiza- 
tion of the grand potential i7 = tr [p {H — p)] (p is a chemical potential used to eliminate 
the particle number constraint Ne = tr [p]) with respect to p amount to an 0{N^) operation. 
In the latter case this a result of the idempotency constraint. 

In order to achieve an 0{N) solution to the problem, LNV use the following strategy. 
First, they take advantage of the fact that the density matrix is local in real spaceBiZi (in the 
sense that pij -^ as Rij -^ oo, where Rij is the distance between basis orbitals 0j and (pj), 
and introduce a trial density matrix X which is set to zero for Rij > R^ {Re is chosen large 
enough to get a good approximation to the true density matrix). Second, no idempotency 
constraint is explicitly imposed; rather, they make use of the purification transformation 
proposed by McWeeneyS: 

p = 3X2 - 2^3 . (3) 

This transformation is such that a matrix which is nearly idempotent {Xx = ^ + S or S, 
\6\ << 1, where Ax is an eigenvalue of X) transforms into a matrix which is more nearly 
idempotent [Ap = 1 — 0{6'^) or +0{6'^)]. Then, p is treated as a physical density matrix 
(i.e., tr[pA] gives the physical expectation value of operator A) and X as a trial density 
matrix whose elements constitute the variational degrees of freedom to be determined by 
minimization of the grand potential 

n = tT[p{H - p)] = tr [(SX^ - 2X^) {H - p)] (4) 

with respect to X. As shown in Ref.EJ, fi in Eq. (^) has a variational local minimum (i.e., 
Q > inexact) at which p is idempotent to second order. Since the number of degrees of freedom 
per atom is fixed by Re and no diagonalization or orthonormalization step is performed, the 
above procedure amounts to an 0{N) solution to the problem. 



We now extend the LNV energy functional to a non-orthogonal basis {0j}, with the 
overlap matrix given by Sij = {(pi \ (f)j). In what follows we use Xij = {(pi \ f> \ 4>j) for 
the trial density matrix in the non-orthogonal basis, to be consistent with the notation 
introduced above. The eigenstates of H are given by Eq. (|1|) and the coefficients {c„j} are 
determined by solving the secular equation 

Y^ {Hij - enSij) Cnj = , (5) 

j 
where Hij = {(pi \ H \ (pj). 

Let C be the matrix defined by Cj„ = c„j (i.e., C has the eigenvectors {V'n} as its 
columns); it then follows that C defines a congruence transformation that diagonalizes H, 
S and X simultaneously: 

C^HC = A , 

C^SC = I , 

C^XC = Xjj , (6) 

where / is the identity matrix, and A^n = ^nSmn and (Xj^)^,^ = 6{ji—en)5mn are, respectively, 
the matrices of H and p in the basis {ipn} of the eigenvectors of H [9{x) is the theta function] . 
From Eq. (^ we have 

5-1 = CC^ . (7) 

Using Eqs. @ and (^, and p^ = 3Xjj — 2Xfj following Eq. (H), we can now generalize 
Eq. (D as follows: 

n= tr [ph (A - /i)] 



tr 



(3S-^XS-^XS-^-2S-^XS-^XS-^XS-^)h'] , (8) 



where H' = H — pS. 

As a matter of convenience, we would like to eliminate S~^ from the energy expression 
in favor of 5*. This can be accomplished by defining the two new quantities 



X = s-^xs-^ , 

p = 3XSX - 2XSXSX , (9) 

as alternative representationsHI for the trial and physical density matrices, respectively. We 
observe that X is a more natural representation of the density operator, in the sense that 
Eq. (|) still holds, i.e., Xij = Y^nCnfinr whereas X,j = EnEkiSikC^kCniSij- Furthermore, 
the expectation value of any operator is given by (A) = tr [XA], where Aij = {(pi \ A \ (pj). 
In terms of X and p the particle number becomes 

Ne = tr [pS] = tr [(3XSX - 2XSXSX) s] , (10) 

and the energy functional is written 

n = tr [pH'] = tr [{SXSX - 2XSXSX) H'] . (11) 



To show that Eq. (11) is variational with respect to X, we note the following. At the 
solution, the density matrix must obey the idempotency constraint Xjj = X^j, which in the 
new representation is expressed as 

XSX = X ; (12) 

furthermore, it must also commute with the Hamiltonian, i.e, Xj^A = AX^j. In terms of X, 
we have 

SXH = HXS . (13) 



Eqs. (|T2D and (|l^) can easily be obtained by applying the transformation generated by C to 
the more familiar expressions in the basis {ip}, and then using X = SXS. From Eqs. (|12D 
and (|I3|) it follows immediately that the variational gradient 
AO 

-J = 3 [SXH' + H'XS) - 2 i^SXSXH' + SXH'XS + H'XSXS) (14) 

vanishes at the solution, thus showing that fl is variational with respect to X. 



Eqs. (p!OD, (|TTD, and (^4|) constitute the central results of this work. Note that the 
standard LNV scheme is recovered from these equations upon substituting Sij = 6ij. 

Before proceeding further, let us comment on the real-space localization properties of 5* 
and 5*"^. We are interested in the case where the basis orbitals are localized in real space, and 
therefore S is also localized. It can be shown that S~^ is then exponentially localizedliSHU, 
with a decay length that depends on the spread of the eigenvalues of S. If S is an ill- 
conditioned matrix, i.e., max(A5)/min(A5) ^ 1 [max(As') and min(As') are, respectively, 
the maximum and minimum eigenvalues of S], then S~^ has a long decay length. 

The advantage of using Eq. (pH]) is that it eliminates the need to compute S^^. A possible 
concern, in making use of the current approach, may be that the matrix X may decay more 
slowly with distance than would the density matrix expressed in terms of orthogonal basis 
orbitals. This may happen when S is ill-conditioned, such that S^^ has a longer exponential 
decay than XO. When this is the case, it becomes necessary to increase the cutoff radius Re 
to obtain the same level of accuracy. However, as we discuss below, our numerical evidence 
suggests that X and X will, in general, have very similar decay as each other, and very 
similar to that of X for an orthogonal basis, and therefore the transformation leading from 
Eq. (^) to Eq. ([TT]) typically preserves the localization properties of Q. Note also that the 
presence of S in Eq. ([TTl) [as compared to Eq. @)] does not affect the linear scaling of the 
method, since 5* is as localized as if in a local basis. 

Next, we present some numerical tests for a three-dimensional tight-binding (TB) model 
for silicon, to illustrate the localization properties of X and X. We use the universal TB 
model proposed by Harrison in its extension to a non-orthogonal basial^. The matrix ele- 
ments of S" in a minimal sp^-basis are given by Su/n = 2kViiin{ti + e^)^^, where / and /' run 
over s and p orbitals and n indicates the type (a or tt) of interaction. The VJ^^'s are the 
universal TB parameters introduced by Harrisonli^, ei are atomic on-site energies and k is an 
adjustable parameter. For the Hamiltonian matrix we have if;;'„ = (1 + A; — ^DVJj'n, where 
S2 is the overlap between two sp^ hybrids ^2 = {Sssa — '^V^Sgpa — ^Sppa)/4:. Both H and S 
are restricted to first neighbors only. For simplicity we set A; = 1 which is very close to the 



value 1.042 commonly used for silicon. 

In Fig. |1] we show the behavior of X and X as functions of the distance i?jj, between two 

j- 2 -I 1/2 

Hap \^ap\ /-^ ' where a and f3 run 



orbitals 0j and (j)j. Plotted is the norm \\X''-^ 



over {s,px,Py,Pz} and A^ = 4 is the block dimension. Also shown is the behavior of Xortho 
for orthogonal orbitals, which is obtained by setting A; = in the TB model. It can be seen 
that within a distance Rij = 8.00 A, both X and X as well as Xortho decay to ~ 2.0% of 
their values at the origin. 

We have also calculated X and X for a basis with an ill-conditioned (almost singular) S, 
by setting k =1.35 (this implies a large overlap between neighboring orbitals). The results 
are shown in Fig. ^ Also shown is the long range behaviour of S^^. Although S~^ has a 
very long decay length in this case, X and X still decay very similarly as Xortho- The point 
here is to show that, provided that the basis is sufficiently localized, X can be as localized 
as Xortho for an orthonormal basis, even when S~^ is ill-conditioned. In any case, by using 
Eq. ([Tl|), the convergence of both quantities, p and S"^, is built into a single quantity X, 
and is controlled by a single parameter R^. 

Because of the fact that S is of the same range as H (i.e., R^ = Rh), the GDM scheme 
preserves the 0{N) scaling of the original method. Nevertheless, the presence of S in 
Eq. (|l^) implies an increase in the scaling pre-factor. In the orthonormal case, the time- 
dominant step involves the calculation of a product Xjk{XH)ki of two matrices of range R^ 
and Re + Rfj-, out to a radius Rij < Re- In the non-orthonormal case, the corresponding 
dominant operation involves two matrices {SX)jk{SXH)ki of range Rc + Rh and Rc + 2Rfj, 
calculated also up to the radius Re [see Eq. (p!4D]. In order to estimate the slowdown factor, 
we determined the ratio of the number of terms that contribute in each case, which was 
found to be 3.6 using the Re and Rj^ of Ref.lil. 

In summary, we have presented a generalization of the LNV density-matrix approach to 
the case of a non-orthogonal basis. An alternative real-space representation of the density 
operator is introduced, which is argued to have similar localization properties as the conven- 
tional density matrix, as suggested by the numerical evidence presented. In this generalized 



energy functional, only the overlap matrix S appears explicitly (as opposed to S~^). The 
new functional is shown to retain its variational property and its linear scaling with system 
size. 
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FIGURES 
FIG. 1. X, X and Xortho SlS functions of the distance Rij between two basis orbitals, for the 
tight-binding model of the text. 

FIG. 2. Same as before for ih-conditioned S {k = 1.35); long range behaviour of S^^ is also 
shown. 
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